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We calculate the conductivity of arbitrarily stacked multilayer graphene sheets within a relax- 
ation time approximation, considering both short-range and long-range impurities. We theoretically 
investigate the feasibility of identifying the stacking order of these multilayers using transport mea- 
surements. For relatively clean samples, the conductivities of the various stacking configurations 
depend on the carrier density as a power-law for over two decades. This dependence arises from a 
low density decomposition of the multilayer band structure into a sum of chiral Hamiltonians. For 
dirty samples, the simple power-law relationship no longer holds. Nonetheless, identification of the 
number of layers and stacking sequence is still possible by careful comparison of experimental data 
to the results presented here. 

PACS numbers: 72.80.Vp,73.23.-b,72.80.Ng 



I. INTRODUCTION 

Enormous progress has been made in the five years 
since the first experiments demonstrated that charge car- 
riers in single graphene sheets behave like massless Dirac 
fermions. Making and studying such single-atom thick 
carbon sheets is now routinely done using a wide variety 
of techniques. (For reviews, see Refs. [HI^. 

While the allure of manipulating single monoatomic 
sheets has understandably attracted most of the atten- 
tion in this field, from a technological, or for that mat- 
ter, fundamental point of view, the properties of few- 
layer-graphene sheets are equally attractive. Many of 
the methods used to make monolayer graphene, such as 
mechanical exfoliation of graphite, epitaxial growth from 
silicon carbide and chemical vapor deposition on metals, 
can be suitably adapted to make graphene stacks, with 
a controllable number of atomic layers. Many of the un- 
usual properties of the Dirac Hamiltonian that are used 
to describe monolayer graphene, such as having chiral 
and ambipolar carriers, survive in these multilayers, so 
one might expect that these sheets would have high mo- 
bility and favorable carrier transport properties. As a 
result, multilayer graphene may play an important role 
in future electronic devices where its additional "layer" 
degree of freedom could be manipulated^ to achieve de- 
sirable properties, such as the demonstrated gate-tunable 
band-gap in graphene bilayersi^ 

There has been experimental and theoretical work on 
the optical properties of graphene multilayers, as well as 
some very recent theoretical predictions on the phonon- 
scattering in these multilayers.^ However, there has not 
been a systematic study of the low temperature transport 
properties of graphene multilayers. In anticipation of 
forthcoming experiments on these systems, we use both 
analytical and numerical methods to understand carrier 
transport in graphene multilayers. 

The complication in studying multilayers is the cou- 



pling between the layers. The carrier transport in a 
single graphene sheet can be readily understood using 
the Dirac Hamiltonian, which is the low energy effec- 
tive theory for vr-orbitals located on the vertices of a car- 
bon honeycomb lattice. For multilayers, the additional 
coupling between orbitals on neighboring layers depends 
sensitively on many factors such as the distance between 
the layers and their relative orientation. For example, 
graphene bilayers with a twist angle between their re- 
spective primitive cells are predicted to largely act as 
decoupled sheetsj^ii For Bernal stacking (also called A-B 
stacking), on the other hand, half the carbon atoms in 
each hexagon of the top layer lie exactly over the center 
of a hexagon of the layer below it. The resulting strong 
coupling between the two layers gives a low energy effec- 
tive theory with a zero-gap hyperbolic dispersion. 

While height fluctuations (or equivalently, having spa- 
tial fluctuations in the interlayer coupling strength), or 
allowing for arbitrary rotations and slips between the 
layers are important for some systems (such as epitax- 
ial graphene), their effects are beyond the scope of the 
present worki^ Here we consider multilayers that come 
in families where the orientation of the upper layers is 
determined by symmetry considerations from the orien- 
tation of the bottom layer. This would be the case, for 
example, if the multilayer inherits its structure from a 
parent structure, as is the case in the mechanical exfoli- 
ation technique. 

We further restrict our multilayer analysis to the lower 
energy stacking sequences in which neighboring layers 
share only one sublattice. For example, we consider bi- 
layer graphene that is A-B stacked, where the two lay- 
ers share one sublattice, but not A-A stacked, where 
each carbon atom of the top layer lies exactly on top 
of a carbon atom of the bottom layer (sharing both sub- 
lattices). The consecutive A-A stacking is energetically 
unfavorable,— so we do not consider this stacking and its 
generalization in multilayer stacks. For trilayer graphene. 



2 



(a) ^^^^^ (b) 




FIG. 1: (Color online) (a) Schematic illustration of (a) three 
types of stacking arrangements, labeled by A, B and C. The 
honeycomb lattice of a single sheet has two triangular sublat- 
tices, labeled by a and /3. (b) Each added layer cycles around 
this stacking triangle in either the right-handed or the left- 
handed sense. Reversals of the sense of this rotation tend to 
increase the number of low-energy pseudospin doublets. 



we consider two possibilities: A-B-A stacking (also called 
Bernal-like) is a Bernal bilayer, with the third layer hav- 
ing carbon atoms located directly above the bottom layer; 
and A-B-C stacking. 

Figure [T] illustrates the different stacking sequences 
graphically. Since the honeycomb lattice of a single 
graphene sheet comprises two interpenetrating triangu- 
lar sublattices, we label the sublattices of each layer a 
and /3. When a subsequent graphene layer is placed on 
top of the stack, we consider the stacking orders where 
either the atoms of the a or the (3 sublattices are dis- 
placed along the edges of the honeycomb of this top 
sheet. This gives a stacking rule that implies three dis- 
tinct but equivalent projections (labeled A, B, and C) of 
the three-dimensional structure's honeycomb-lattice lay- 
ers onto the x-y plane and consequently 2^~^ distinct 
stacking sequences for an A^-layer stack. 

The electronic properties of multilayer graphene 
strongly depend on the stacking sequence. Periodically 
stacked multilayer grapheneiiiiS and arbitrarily stacked 
multilayer grapheneii have been studied theoretically, 
demonstrating that the low-energy band structure of a 
graphene multilayer consists of a set of independent pseu- 
dospin doublets. It was shown that an energy gap can 
be induced by a perpendicular external electric field in 
ABC-stacked multilayer graphenei ^'^'^^ Furthermore, in 
ABC stacking, electron-electron interactions play a more 
important role than other stacking sequences due to the 
appearance of relatively flat bands near the Fermi level 
This enhanced role of electron interactions raises the like- 
lihood of strongly-correlated ground-states, a possibil- 
ity that we ignore in our semiclassical treatment below. 
Optical properties of multilayer graphene using absorp- 
tion spectroscopy have been studied experimentally^^ 
and theoreticall y-^ ''i-^^ showing characteristic peak posi- 
tions in optical conductivity depending on the stacking 
sequence. 

Transport properties of monolayer, bilayer and mul- 
tilayer graphene have been studied theoreticallyi^Ti^^ 



within coherent potential approximations. These ap- 
proaches capture the scattering properties of the impu- 
rity potential (which is important for strong disorder), 
but they are often restricted to small system sizes, and 
do not accurately account for the disorder-induced spa- 
tial inhomogeneity of the fluctuating local carrier density. 
We believe this inhomogeneity dominates the transport 
properties of these graphene multilayers (see discussion 
in Ref. 2). 

Our main finding is that for relatively clean samples, 
the carrier density dependence of the multilayer conduc- 
tivity follows a power-law dependence for more than two 
decades, a direct consequence of the effective low energy 
chiral decomposition. For dirty samples, the carrier den- 
sity inhomogeneity induced by the disorder washes away 
this power-law relationship. However, the various stack- 
ing sequences give characteristically different dependence 
of the multilayer conductivity on carrier density. By care- 
ful comparison with experimental data, our results could 
be used to identify both the number of layers and the 
stacking sequence of a multilayer graphene sample. 

The rest of this manuscript is organized as follows. In 
Sec. ini we describe the theoretical model where we solve 
for the multilayer graphene band structure using a tight- 
binding model that includes both the nearest-neighbor 
intralayer hopping and the nearest-neighbor interlayer 
hopping, and solving for the conductivity within the 
Boltzmann transport formalism. In Sec. IIIII we present 
our results for graphene stacks comprising one, two, three 
and four layers, treating impurity scattering by both 
Coulomb potentials and short-range disorder. In the ap- 
pendices we present details of the chiral decomposition 
and the transport properties of J-chiral fermions, as well 
as analytic results for the electronic transport in bilayer 
graphene. 

II. THEORETICAL MODEL 

A. Tight binding Hamiltonian 

The low energy effective Hamiltonian for the tt- 
orbital continuum model for arbitrarily stacked A^-layer 
graphene centered at the hexagonal corners of the Bril- 
louin zone is given byS2ii^ 

n = J2^lHip)^p, (1) 
p 

where = {ci,a,p, ci,/3,p, • • • , CN,a,p, CN,p,p) and c;,^,p is 
an electron annihilation operator for layer Z = 1, • • • ,N, 
sublattice /i = a, /3 and momentum p measured from K 
or K' point. 

The simplest model for a multilayer graphene system 
allows only nearest-neighbor intralayer hopping t and 
the nearest-neighbor interlayer hopping t±_ . The in-plane 
Fermi velocity for monolayer graphene, vq, is related to 
t by ^ = ^t, where a = 0.246 nm is the lattice con- 
stant of monolayer graphene. This model ignores some 
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aspects of the electronic band structure - in principle, 
corrections to this model, such as adding next-nearest- 
neighbor hopping, can easily be included, although in 
practice, it is often numerically quite intensive. We find 
that such corrections do not significantly alter any of our 
main findings. 



B. Boltzmann transport theory 



in Appendix [Cj For convenience, we define the dimen- 
sionless potential Vimp(?!') 

efcp Kp 
where e is the dielectric constant and a 



ehvQ 



is the 



effective fine structure constant. The conductivity can 
then be written as 



The conductivity is a property of electrons close to 
the Fermi energy and is obtained from the Einstein re- 
lation, CTB = e^'D{eF)D where I'(eF) is the density of 
states at the Fermi energy ep and D is the diffusion con- 
stant. In cases in which the Fermi surface has multiple 
sheets (lines in the two dimensional cases we consider 
here), the conductivity is the sum over such contribu- 
tions for each sheet. See Sec. Ill Dl for details. For con- 
venience, we denote the density of states for each sheet 
as 25(ep) = (7s3vP(eF) where gs = 2 and = 2 are spin 
and valley degeneracy factors. For diffusive transport in 
two dimensions, D = ^fpTp. The Fermi velocity, density 
of states and relaxation time can be calculated from the 
dispersion relation as 
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where V^^p is the squared effective impurity scatter- 
ing potential averaged over the angle, and the Fermi 
wavevector is related to the carrier density and applied 
back gate voltage fc| = 47rn/(gs5v) oc Vq- When the 
Fermi surface has multiple sheets, the left hand side of 
this last relation should be summed over kp for each 
sheet. The validity of the Born approximation implicit 
in this formulation is discussed in Appendix |^ and the 
special case of J-chiral fermions is treated in Appendix IB] 
The scattering matrix element V^^^p that gives rise to 
the transport relaxation time is obtained using the Boltz- 
mann transport formalism 
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where Vinip('/') is the matrix element of the impurity po- 
tential at scattering angle 4>, and F{(j)) is the chirality 
factor that arises from the projection of the spinor wave- 
functions between the incoming and outgoing states. For 
the case of intervalley and interband scattering, the treat- 
ment of the chirality factor is more subtle and discussed 
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where n — gggvkp / {An) is the carrier density, and all the 
information about the band structure and the nature of 
the disorder potential is captured by v-p and VJ^p. When 
several bands cross the Fermi energy, we calculate v-p and 
^mp (which yields ai) for each of the i bands and then 
calculate the total conductivity as a = where the 

applied gate voltage is proportional to n = ^i^i, and 
ni — gsgvkp i/(47r) is fixed by keeping ep the same for all 
bands. 



C. Impurity Scattering 

Different types of impurity potentials give qualita- 
tively different results for the conductivity. This can 
be seen in Eq. [5l where the wavevector dependence 
of the Fourier transform of the impurity potential 
Mmp[9 = 2kpsm{(j)/2)] changes the scaling of the con- 
ductivity cr(n). Studies on monolayer graphene have 
explored a wide variety of disorder potentials including 
long-range Coulomb (T4mp(9) 9~^)i Gaussian-white 
noise (V^mp('?) ^ 9°) and Gaussian-correlated disorder 
(V^mp('?) ^ exp[— g^]) as well as resonant scatterers that 
cause a maximal phase shift of tt/2 between incoming 
and outgoing wavefunctions. We refer to Ref. |3 for more 
details on the different kinds of disorder in graphene. 

For our purposes, we focus on what we believe to be 
the most relevant scattering mechanisms, i.e. charged im- 
purities (which act as screened long-range Coulomb dis- 
order) and short-range defects (approximated as white- 
noise disorder). Using the Thomas-Fermi screening the- 
ory, one can write expressions for the dimensionless scat- 
tering potential which was defined in Eq. U) For 



screened charged impurities from V^mp(0) 
with q = 2kp sin((/)/2) we find 



-qdh 



imp 







d(t) F(0)(1-COS0) ^_4fcFd.„„ sin(0/2) 

27r (2sin(0/2) -H^tf)^ 



(6) 

where the dimensionless Thomas-Fermi wavevector is 
gxF — iTp/kp = gsgvOt{vQ / vp) and dimp is the distance 
between the impurities and the graphene sheet. When 
several bands cross the Fermi energy, the Thomas-Fermi 
screening wavevector determined from the total density 
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of states including all the bands. The monolayer Fermi 
velocity Vq was defined below Eq. [1] 

The short-range disorder potential can be character- 
ized by an effective scattering cross-section length dgc 
such that Vi„,p{(t)) = ^^^s!^ or 
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Taken together with Eq. [SJ this completely defines the 
electrical conductivity in terms of the multilayer band 
structure discussed in Sec. Ill Al 



D. Intervalley and Interband Contributions 

At typical carrier densities, the band structure of 
monolayer graphene comprises two Dirac cones that are 
centered at two inequivalent points (also called valleys) 
of the Brillouin zone boundary labeled as K and K'. The 
scattering of carriers between the valleys requires a mo- 
mentum transfer of Q = 47r/(3a) ss 17 nm~^, and is 
strongly suppressed for Coulomb impurities. For short- 
range scattering, the treatment depends on how one mod- 
els the intervalley scattering matrix element, but in most 
cases it is sufficient to consider a single valley with a 
suitably defined impurity concentration nimp (see e.g. 
Ref. [251 ). For concreteness, we assume that for short- 
range impurities, the intervalley and intravalley scatter- 
ing matrix elements are equal, that riimp is the aver- 
age impurity concentration in a single valley (and is the 
same for both valleys). The conductivity using these 
definitions is smaller by a factor of 2 from the case of 
n-imp = ntot = UA + nB, where ua and ub are concentra- 
tion of impurities or defects on the A and B sublattices, 
respectively. 

As discussed in Appendix [B] within a single valley, all 
finite stackings that are a subsequence of periodic ABC 
stacking; i.e. A, AB, ABC, ABCA, etc.; have only a sin- 
gle band at low energies, and the chirality increases as 
the number of layers increases. However, for other stack- 
ings, the band structure features multiple chiral bands 
with different dispersion relations that are centered at the 
Dirac point. For Coulomb impurities, the matrix element 
in Eq. [5] can be computed for both intraband and inter- 
band scattering, and their scattering rates added in ac- 
cordance with Matthiessen's rule. For short-range impu- 
rities, the interband scattering is more subtle. A straight- 
forward application of Eq. [7] would lead to a strong sup- 
pression of all interband scattering because the chirality 
factor F{(f)) vanishes or is significantly less than one. As 
discussed in Appendix [Cl we believe this to be unphysical 
because it requires the short-range impurity potential to 
be diagonal in the space of all the layer and valley quan- 
tum numbers. Most short range scatters we can imagine 
would be localized to one layer and a specific sublattice 
so the scattering potential would not be diagonal in this 
space. The screened Coulomb potential, on the other 



hand, varies slowly between the layers and sublattices 
and can be much better approximated as diagonal in 
that space. In the absence of a microscopic model for 
a particular impurity model (e.g. computing the scatter- 
ing potential resulting from a single vacancy or from the 
binding of a single hydrogen atom to the top layer), we 
believe that for short range scatterers, a more realistic 
assumption is to set F{(j)) = 1. This correctly weights 
the relative importance of intraband and interband scat- 
tering, and therefore gives the correct qualitative carrier 
density dependence of the conductivity. 

The role of interband scattering is most striking when 
one considers the large density regime where higher en- 
ergy bands become accessible. For screened charged im- 
purities, the additional density of states in these higher 
energy bands enhance the screening of long-range impu- 
rities which will sharply increase the conductivity, while 
the interband scattering is suppressed by the chirality 
factor. On the other hand, for short-range impurities, the 
higher energy band becomes an additional source of in- 
terband scattering that sharply decreases the conductiv- 
ity. At the time of writing, there have been no transport 
experiments that could probe these higher energy bands 
by inducing sufficiently large carrier densities. However, 
if such an experiment is done in the future (perhaps by 
finding better electrolytes), then the increase (decrease) 
of (j{n) would be indicative of Coulomb (short-range) im- 
purities being the dominant source of scattering. 



E. Effective Medium Theory (EMT) 



At low carrier density, the disorder induced fluctua- 
tions in the local density become larger than the spa- 
tially averaged carrier density. This has been called 
the electron-hole puddle regime. We use the effective 
medium approach developed in Ref.l26lto obtain the bulk 
conductivity ct^mt of this inhomogeneous medium. It was 
shown in Ref. 27 that by assuming a Gaussian probabil- 
ity distribution for the carrier density, (Temt('T') could be 
obtained from o^{n) using 



I dn' exp 
Jo 



o■B(f^') 



cosh — — — — -— = 0, 

(8) 

where nj-ms parameterizes the Gaussian distribution, and 
aB{n) is obtained numerically from Eq. [5l This in- 
tegral equation extrapolates from crEMT(?T-) w min 

for 

I"- 1 ^ n^ms to a^^„:{n) « a^{n) for \n\ > rirms. The 
effect of the puddles, therefore, is to give rise to a min- 
imum conductivity plateau where the conductivity re- 
mains roughly constant when the average of the carrier 
density is smaller than its fluctuationsi^ 
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III. RESULTS AND DISCUSSION 

As outlined above, for an arbitrary graphene stack we 
first solve Eq. [T] to obtain the band structure. For sim- 
plicity, we choose t = 3 eV, t± = 0.3 eV, a = 1 and 
n-ims = JT-imp- Neglecting higher order hopping terms 
keeps the spectrum rotationally symmetric (reducing the 
computational time), while the ratio of nmis/'T-imp is a 
number of order unity that can be calculated within the 
low-density chiral decomposition (see Appendix IB|) . We 
do not believe that these approximations significantly al- 
ter our findings. 

Having solved numerically for the wavefunctions, one 
can compute vp fEg. Bal) . p{e) fEa. l2b|) . as well as the chi- 
rality factor (see Appendix [Ct. For the short-range 
and Coulomb disorder potential we use representative^^ 
parameters: dsc = 0.3 nm and dimp — 1 nm. Equation [5] 
then gives the predicted carrier density dependence of 
the conductivity. 

Figure [5] shows the results for short-range scatterers 
such as point defects. The left panel (Fig. [5^ and Fig.[5]D) 
assumes a relatively clean sample with njmp — 10^*^ cm~^. 
The solid lines are CTemt calculated from Eq. |8l while 
dashed lines show the (approximate) power-law depen- 
dence of the conductivity on carrier density. We find 
that for most cases the conductivity limited by short- 
range scatterers exhibits a unique power-law (for more 
than two decades) that depends on the number of lay- 
ers and the stacking sequence. For such clean samples, 
where the scattering is dominated by short-range disor- 
der, with the exception of the similarity between AB bi- 
layers and ABA trilayers, the distinct power-law depen- 
dence for a{n) enables easy identification of the sample 
from transport measurements. 

The right panel (Fig. [2j; and Fig. [2ji) shows the re- 
sults for a dirty sample (riimp = 10^^ cm~^). The 
solid lines show the EMT result (Eq. [5]) and the dashed 
lines show the Boltzmann result before the EMT av- 
eraging (Eq. [5]). We note that the results for mono- 
layer and bilayer graphene agree with analytical calcu- 
lations discussed in Appendix [Dl While there is no sim- 
ple power-law dependence of the conductivity on carrier 
density, the different stacking sequences have very differ- 
ent functional forms for a{n). (For example, the ratio 
a{n = 10^^ cm~^)/(T(n — 10^^ cm~^) varies by almost 
an order of magnitude). It might therefore still be pos- 
sible to identify the number of layers and the stacking 
sequence from transport. 

The results for charged impurities is shown in Fig. |31 
Again the left panel shows a clean sample (nimp = 
10^" cm~^). While all the different stacking sequences 
have a power-law dependence for cr(n), unfortunately, 
with the exception of ABA and ABCB, the rest are 
all very close to a linear dependence. It would there- 
fore be difficult to distinguish the samples in any trans- 
port measurement in the low carrier density regime (i.e. 
n < 10^^ cm~^). For a dirty sample with a larger range 
of carrier density (nimp = 10^^ cm~^,n < 10^'^ cm~^). 



we find that the functional forms for cr(n) are sufficiently 
different. We emphasize that the figures show the con- 
ductivity on a log-log scale. Different slopes correspond 
to different scaling exponents 7, where cr(n) ^ rf . For 
example, as seen in Fig. [3jl, while the different stackings 
of the tetralayer have similar values for the conductiv- 
ity at n = 10^'^ cm~^, the minimum conductivity for the 
ABCB tetralayer is two orders of magnitude larger than 
the ABCA stacking. 

To further understand these results, we note that the 
transport properties of graphene multilayers are deter- 
mined by two characteristic densities: the band density 
riQ = Qsg^k^ I (At:) where hvoko = t±/2, and the impu- 
rity density n-i^p. By gating graphene, one can change 
both the carrier density and the type of carriers, where a 
negative back-gate voltage induces holes, and a positive 
back-gate induces electrons. With special dielectrics one 
can induce carrier densities as large as n « 10^* cm~^ 
(see Ref. HO) , although typically carrier densities do not 
exceed n « 5 x 10^^ cm^^. 

For carrier densities much lower than the characteristic 
band density no ~ 2 x 10^^ cm^^, one can decompose the 
electronic structure of an arbitrary multilayer into par- 
allel channels of bands, each with the simple dispersion 
relation Ck ~ k'^ where J is the chirality index. The num- 
ber of channels and the chirality index are determined 
from the stacking sequence as discussed in Ref. [3l|. More 
details on the wave-functions and transport of J-chiral 
fermions can be found in Appendix |B] 

In the opposite limit of n 3> no, the band structure of 
iV-layer graphene decomposes into that of N decoupled 
monolayer graphene sheets, irrespective of the stacking 
sequence. Since the transport of monolayer graphene has 
been well studied (see Ref. 0), we do not explore this 
limit in any detail. Although we note, that even at the 
extremely large carrier density n w 10^^ cm~^, one is not 
yet in the limit of essentially decoupled sheets. 

The second important scale is that of disorder. Typ- 
ical values of njmp vary from 10^*^ cm~^ (in suspended 
graphene) to 5 x 10^^ cm~^. Only when the carrier den- 
sity is much larger than riimp can one use the usual semi- 
classical Boltzmann transport theory. When n < nimp, 
the inhomogeneous landscape of electron and hole pud- 
dles gives rise to a saturation in the conductivity ap- 
proaching a finite (Tmin at the Dirac point. 

One can immediately identify two regimes that are ex- 
perimentally relevant. When ni,„p <^ n <^ uq, one can 
exploit the chiral decomposition to obtain a power-law 
dependence of the conductivity on carrier density. This 
is what we called the "clean sample" regime (the left 
panels of Fig. [5] and Fig. [3]) . In Appendix |B] we derive 
the Effective Medium conductivity for arbitrary J-chiral 
fermions, and we can use this decomposition to under- 
stand our results. 

As an example, consider the clean ABC stacked tri- 
layer. As seen in Fig. [51 c('t-) ~ for short-range 
impurities. This dependence follows directly from the 
low energy chiral decomposition discussed in Appendix [B] 
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FIG. 2: (Color online) Graphene conductivity assuming short-range scatterers. For clean samples, (nimp 10^'' cm~^), the 
conductivity a-s,mT(ri) has about two decades of power-law dependence, (a) The conductivity (solid lines) for monolayer, bilayer 
and ABC trilayer graphene each follow different power-laws (dashed lines), (b) The conductivity (solid lines) for the different 
stacking sequences for tetralayer graphene also have different power-laws (dashed lines). This indicates that for short-range 
impurities, in most cases transport measurements can be used to identify the number of graphene layers. For dirty samples (c) 
and (d), where riimp ~ 10^^ cm~^, the density dependence is no longer given by a power-law. Solid lines show the conductivity 
after using an effective medium theory to average over the disorder induced carrier density fluctuations. Dashed lines are the 
results before such averaging. Although the conductivity does not have power-law dependence on carrier density, the transport 
properties still strongly depend on the number of layers and their relative stacking-order. Therefore, transport measurements 
could still be used to identify the type of graphene multilayer. 



where the ABC trilayer is approximated by a J = 3 
chiral system, and aj ^ n"'^^ (see Eq. [B5}. Similarly, 
the numerical results for the ABC stacked trilayer with 
charged impurities shown in Fig. [3] is quite close to the 
expected power-law a[n) ^ n. This small difference be- 
tween the numerical results and those anticipated from 
the chiral decomposition is the result of our using a fi- 
nite dimp = 1 nm for the distance of the Coulomb im- 
purities from the graphene sheet. The finite dimp softens 
the small-distance divergence of the Coulomb potential, 
thereby increasing the conductivity and giving a larger 
coefficient 7 for the (approximate) power-law a ^ n'^ . 

A similar analysis can be done for the ABA stacked 
trilayer graphene. At low energies, it is described by 



a direct product of J = 1 and J ~ 2 chiral systems. 
At low density, the band structure has two sheets, one 
with a linear dispersion like monolayer graphene, and one 
with a parabolic dispersion similar to bilayer graphene. 
Requiring a constant Fermi energy, one can introduce 
a dimensionless parameter x = ep/t^. The chiral de- 
composition is valid when x ^ 1. One notes that 
X {v* / vqY {k-p / ko)'^ where the effective Fermi velocity 
V* ~ vq is shown in Table ID This implies that at low car- 
rier density where the chiral decomposition is valid, the 
band with the larger chirality also has a larger carrier 
density and larger density of states. In the case of ABA 
trilayers, this implies that for short-range scatterers, it 
behaves exactly like the AB bilayer which is also a J = 2 
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FIG. 3: (Color online) Graphene conductivity (solid lines) assuming charged impurity scattering, (a) For clean samples, the 
long-range Coulomb scatterers give similar power-law dependence (dashed lines) of the conductivity on carrier density for 
monolayer, bilayer and ABC trilayer graphene. (b) The ABCA and ABAB stacking orders of tetralayer graphene give a similar 
power-law dependence (dashed lines) to monolayer and bilayer graphene. We conclude that the low density transport properties 
for many graphene multilayers look quite similar under long-range scattering, making it difficult to distinguish between them 
in a transport measurement, (c) For dirty samples, the transfer curves for monolayer, bilayer and trilayer graphene are each 
quite different and can be easily distinguished. (Notice the logarithmic scale on the y-axis). (d) For dirty tetralayer graphene, 
the conductivity depends strongly on the stacking sequence, and transport measurements could distinguish between the various 
types of stacking. 



chiral system and crj=2 ~ n (as seen in Fig. [5^). For 
charged impurities, the large density of states from the 
J — 2 band effectively screens the impurities so that the 
J = 1 band behaves like monolayer graphene with short- 
range impurities having a ^ constant. All the shown 
power laws in left panels of Fig. [2] and Fig. [3] can be un- 
derstood in this manner using the chiral decomposition. 

The second regime relevant to experiments is when 
'^imp ^ «o ^ max(n). We called this the "dirty sample" 
regime, since having a cleaner sample offers no qualitative 
difference. The important point is that since n > uq, the 
chiral decomposition is not valid, and the band-structure 
has no simple analytical form. This regime can be seen 
in the right panels of Fig. [5] and Fig. [3] Although no 
simple analytical expression or power-law behavior deter- 



mines cr(n), the results for the different number of layers, 
the different stacking orders and short-range vs. long- 
range impurities are all sufficiently different. Therefore, 
by comparison of experimental data to the results pre- 
sented here, it should be possible to identify not only 
the number of layers and stacking sequence, but also the 
nature of the dominant source of disorder. 



IV. CONCLUSION 

We have considered the transport properties of mul- 
tilayer graphene stacks. The formalism is quite general 
and can be used for TV-layers of graphene with arbitrary 
stacking between the layers. In the absence of any exper- 
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imental data for layers with A'' > 2, we have considered 
the most energetically favorable stacking sequences and 
the cases of both short-range and long-range impurities. 
For monolayer and bilayer graphene, our results agree 
with previously known results (see Appendix |D]). For 
trilayer graphene, we show that ABA and ABC stack- 
ing have very different transport properties and can be 
distinguished from each other. Similarly, for tetralayer 
graphene, ABCA, ABCB and ABAB each is a sepa- 
rate electronic material with its own characteristic de- 
pendence of conductivity on carrier density. (The ABCB 
and ABAC stackings have the same conductivity since 
they are related to each other by a symmetry transfor- 
mation.) 

An important objective of this work is to enable exper- 
imentalists working on multilayer graphene to be able 
to use transport measurements to identify and charac- 
terize their multilayer samples. In addition, one could 
use our results to determine the nature of the dominant 
impurity potential^ the effect of changing the dielec- 
tric environment^^ or identify when other effects (such 
as quantum interference^) that we have neglected in our 
semi-classical approach become important. 

Information from transport measurements could be 
used in conjunction with other techniques such as Ra- 
man spectroscopy'^*, optical detection^^ and observing 
the scattering from phonons at high temperature."' Our 
main finding is that in most cases, the different stacking 
sequences have different electronic properties that result 
in characteristic dependence of the conductivity on the 
applied carrier density. 
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impurity potential determines the conductivity. For con- 
sistency with the notation in Ref. 35, we introduce a scat- 
tering potential Mmp = Vo in Eq. [2c] so that Eq. [5] reads 



9s9v 



1 



2 h 



(A2) 



In the Born limit, it is the product TiimpFo^ that deter- 
mines the conductivity. 

Ultimately, microscopic models of measured defects 
will determine whether the scattering is better described 
by the unitary limit or the Born limit. In the absence 
of such a model, an important issue for our results is 
whether the combination of conductivities and carrier 
densities we consider require that the scattering poten- 
tial be so strong that the Born approximation is no longer 
valid. In a recent article, Ferreira et al^ test the valid- 
ity of the Born approximation by computing the scat- 
tering amplitude as a function of the scattering potential 
strength for short-range impurity scattering in monolayer 
and bilayer graphene. For weak scattering potentials, 
the Born approximation results agree with their more 
general calculations but for large enough potentials, the 
scattering amplitude reaches the unitary limit where the 
scattering phase-shift is 7r/2 irrespective of the strength 
of the impurity potential. We use their results to argue 
that the conductivities and charge densities we treat are 
consistent with the Born approximation being valid. 

As discussed by Ferreira et al. the validity of the 
Born approximation depends on the quantity A = 
{yQ/'KhvY)kY\n.{k-pR). If ^ ^ 1, we have the unitary 
limit, and if A ^ 1 we are in the Born limit. To check 
the self-consistency of the Born approximation limit, we 
can re-write this condition as 
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Appendix A: Validity of the Born approximation 

In this work, we use the Born approximation for both 
short-range and Coulomb impurities as seen explicitly in 
Eq.[7l This approximation is valid if the scattering poten- 
tial is weak enough. As the scattering potential becomes 
stronger, the scattering approaches the unitary limit, in 
which the scattering amplitude becomes independent of 
the strength of the potential. In the unitary limit, the 
strength of the impurity potential drops out of the ex- 
pression for the conductivity, which then depends only 
on the impurity concentration^ 



4e2 



(J = 



h 2Trh 



■\n{kFR). 



(Al) 



'imp 



where R characterizes the range of the potential)^ In 
the opposite limit of weak scatterers, the strength of the 



A' = 



Vn 



TThv\ 



■fcp In(fcF-R) 



4 n e^lh 
2 '— < 1, 

TT nimp cr 



(A3) 



where we assume ln(A:F^) ~ 1. A similar expression is 
obtained for bilayer graphene with the prefactor S/tt re- 
placed by 1/ (47r) . The Born approximation is valid for a 
relatively high density of relatively weak scatterers. Since 
a > Ae^ /h (both experimentally, and for the validity of 
the diffusive approximation, see Ref. Issl ) , and typical car- 
rier densities range from 10^° cm~^ to 5 x 10^^ cm~^, the 
Born approximation provides a consistent solution when 
there are no fewer than one short-range impurity per 
40 nm^ (or more than one defect per 2000 carbon atoms) . 
These numbers seem quite reasonable, given the prepa- 
ration method of these samples. In Ref. [2^ the authors 
measured a — 280 /h for monolayer graphene, giving 
two orders of magnitude wider range of carrier densities 
for the validity of the Born approximation for the same 
impurity concentration. Similarly, in Fig.[2^ although we 
take riimp = 10^° cm~^ in order to illustrate the power 
law dependence of the conductivity due to the chiral de- 
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TABLE I: J-chiral decomposition for monolayer, bilayer, trilayer and tetralayer graphene with different stacking arrangements 
(see Ref. [sj for more details). 



Number of layers (N) 


Stacking sequence 


Chirality (J) 


Effective velocity v* /vq 


1 


A 


1 




1 


2 


AB 


2 




1 


3 


ABA 


102 




10 2-1/^ 


3 


ABC 


3 




1 


4 


ABAB 


202 




/(^+1) /2 1/^^(^-1) /2 


4 


ABCA 


4 




1 


4 


ABCB 


103 




l0\/2/2 


4 


ABAC 


103 




l0\/2/2 



composition, the Born approximation is still valid be- 
cause we also take a/{e'^/h) to be much larger. While 
such arguments do not rule out the possibility that the 
scatterers are in the unitary limit, they nonetheless estab- 
lish the Born approximation treatment is self-consistent. 



graphene and J = 2 bilayer graphene, and in general for 
periodic ABC stacking). 

Throughout the manuscript we have used the following 
properties of the energy levels and wavefunctions for the 
J-chiral system 



Appendix B: Conductivity of J-chiral Fermions 

At very low carrier density, an arbitrarily stacked 
graphene multilayer can be described as a superposition 
of pseudospin doublets. This decomposition holds so long 
as hvok-p <C t_L, where J is the chirality index for the 
pseudospin doublet. 

The rules for the decomposition are as follows: (mono- 
layer graphene) A— > (J = 1); (bilayer graphene) AB— >• 
(J = 2); (trilayer graphene) ABA-> (J = 1) ® (J = 2) 
and ABC — >■ (J = 3). This notation indicates, for ex- 
ample, that an ABC stacked trilayer is described by a 3- 
chiral Hamiltonian, while an ABA stacked trilayer is com- 
posed of two bands - one similar to monolayer graphene 
and the second similar to bilayer graphene (see Table U 
and Refs. IsTIIstIIssI for more details). The J-chiral Hamil- 
tonian is defined as 



n = tA 











(Bl) 



where ly^ = hv*ke^'^'' /t± and v* is the effective in-plane 
Fermi velocity (for example, v* — vq for J = 1 monolayer 



£s,k 



St, 



hv*k 



(B2) 



V2 V e 



The band- index s = ±1 corresponds to the positive (neg- 
ative) energy states of the conduction (valence) band and 
ep — £s.k=k-p ■ The intraband chirality factor is then cal- 
culated as 

F(0) = fc, </, = 0|s, fc, 0) |2 = 1 [1 + cos( J0)] . (B3) 

The scaling of the conductivity with carrier density can 
be immediately obtained by noticing that vy ~ A:p~^ and 
p(e_F) ^ fcp^"'. This gives 



,.7-1 



"imp Vimp 



(B4) 



which depends on the scattering potential Vlmp- Assum- 
ing dimp = for simplicity and restoring the dimensions, 
we find 



— { " ^ 



( hv* \ Aim 



Airn 



J-2 



1 



7ra 7,7 



16 

^I3j 



Short-range disorder. 
Bare Coulomb, 

Overscreened Coulomb (a ^ 1), 



(B5) 



where /3j = 1/2 for J = 1, /3j = 1 for J > 1 and /3j = 2 for F(0) = 1, while 7,7 = 1 for F((/)) in Eq. and 
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7j = 2 for F{(j)) = 1. (Here we are considering a J- 
chiral system in Eq. IBll and did not include intervalley 
scatterings.) The result for the bare Coulomb potential 
was presented in Eq. IB5I for a pedagogical reasons, and 
F{(j)) — 1 case for the bare and screened Coulomb poten- 
tials was also considered for completeness. We note that 
since 



we find that for the J-chiral Hamiltonian (with J > 2) 



Qtf 



4a 

T 



(B6) 



for the low density limit hp ^ 0, the overscreened 
Coulomb potential becomes a good approximation for 
charged impurities. In this low density limit, among the 
short-range scattering and screened Coulomb scattering, 
screened Coulomb scattering dominates over short-range 
scattering for J < 2, (with a corresponding a{n) ~ n); 
while short-range scatterers dominate for J > 2 (and 



cr(n) 



^). (Note that because of the Matthiessen's 



rule, a scattering mechanism with smaller conductivity 
dominates.) Bilayer graphene at low carrier density (or 
J = 2) is unique in that both charged impurities and 
short-range disorder give conductivities with the same 
carrier density dependences^— making them difficult to 
distinguish experimentallyjS In the opposite limit of very 
high carrier density, the energy band structure of multi- 
layer graphene separates into a collection of decoupled 
monolayer graphene bands. 1^ As a result, the conduc- 
tivity scales as that of a monolayer at very high carrier 
density. 

As discussed above, the chiral decomposition works 
only at low carrier density where it is known that density 
fluctuations dominate the transport properties 
One must therefore estimate whether there is a regime of 
validity where the carrier density is large enough so that 
the puddle physics no longer dominates (i.e. n ^ nmis), 
but small enough that the chiral decomposition is still 
valid (i.e. n ^ no); here rij-^s is the root-mean-square 
fluctuation in carrier density induced by the disorder po- 
tential, while no ~ 2 X lO^'^ cm^^ is the crossover density 
scale. We estimate nrms using the self-consistent approx- 
imation of Refs. 28! and where (el) = «-imp(^D), and 
Vd is the disorder potential of screened charged impu- 
rities located at some distance dimp from the graphene 
sheet 



{2n) 



2ne^ exp(-g(ii„ip) 



n 2 



e{q + ^tf) 



27rniinp ( ftwo ) ^ C"?^ ( 2gTF rfimp ) 



(B7) 



For the relevant limit n ^ no, Cq {x) 



(for details, 

see Ref.|28|). Although we used a charged impurity model 
for the disorder potential, in this limit the impurities 
are perfectly screened, and the long-range nature of the 
impurity becomes irrelevant (i.e. one gets similar results 
starting from a short-range impurity model) . Assuming a 
Gaussian probability distribution for the carrier density. 



/ Snji^^ipo/^^ 
32^' 



(B8) 



which shows that for any given J, one can determine 
the minimum disorder concentration nimp necessary to 
ensure that n 3> nrms- (The factor of 3 inside the square- 
root is added to conform to the convention for graphene 
monolayers, see Ref. 0.) Note that for J = 1 monolayer 
graphene, one has to use the full dielectric function for 
the screening 1^ 

From Eqs. [8] and IB5I one can easily construct an Ef- 
fective Medium Theory for the J-chiral model. The con- 
ductivity is obtained by solving the integral equation 



/•OO 

/ dn' exp 

^0 



2n2 



cosh 



- CrEMT(") _ Q 
f7j(n') + CrEMT("-) 

(B9) 

To illustrate the differences between different J-chiral 
models, in Fig. 0] we show (7BMT/min((TBMT) for J = 2,3 
and 4 assuming short-range disorder, where we estimate 
nrms from Eq. lB8l assuming that djmp ~ 1 nm and njmp = 
10^^ cm^^. As seen in the figure, the electron and hole 
puddles tends to pin the conductivity value close to its 
minimum value, and that the puddle regime increases 
with increasing J. 



Appendix C: Chirality factor for intervalley and 
interband scattering 

It is often argued that monolayer graphene has a high 
mobility because the chiral nature of carriers prevents 
backscattering. In the diffusive regime, the scattering 
rate involves an integral of the chirality factor F{(f)) over 
all angles weighted by the Boltzmann factor 1 — cos0. 
As discussed in Ref. '39', the enhancement due to chiral- 
ity is no more than a factor of order unity. Similarly 
for graphene multilayers, one can calculate -F(0) numeri- 
cally, and illustrative examples are shown in Fig.[5]for the 
ABA trilayer and ABCB tetralayer. For intraband scat- 
tering and at low density, the chirality factor agrees with 
the analytic results in Eg. IB3I derived for J-chiral Hamil- 
tonians. For the range of carrier densities we consider, 
the chirality factor changes the conductivity by a factor 
of order unity, just like the case of monolayer graphene. 
However, a similar calculation of the interband chirality 
factor shows quite different results. Figure [5] shows that 
the interband chirality factor is exactly zero for the ABA 
trilayer and strongly suppressed for the ABCB tetralayer. 

It can be shown that the interband chirality factor 
in all periodic AB stackings vanishes from the form of 
the wavefunctions^i, and in all other layer stackings it is 
strongly suppressed compared with the intraband chiral 
factor. At first glance, this might suggest that interband 
scattering is negligible and can be ignored. However, we 
point out that being able to decompose the scattering 
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FIG. 4: (Color online) Effective medium theory result (solid 
lines) for the conductivity assuming short-range disorder as 
a function of dimensionless density n/nimp for the J-chiral 
Hamiltonian (see Eq. IB9|I for J = 2, 3 and 4. The results 
assume an impurity density of nimp — 10^^ cm^^. Note that 
o /(Jmin is independent of dsc because dsc appears as a multi- 
plicative factor in Eq. IB5I Dashed lines show minimum con- 
ductivity (i.e. (j/amin = 1), and the Boltzmann conductivity 
(Eq. IB5[) without performing the EMT average. Notice that 
the puddle regime (marked with arrows) gets larger with in- 
creasing J. 



chirality factor F((/)) as we did in Eq.[3]relies on the impu- 
rity potential being diagonal in the 2N x 2N chiral-basis 
for N graphene layers with 2 valleys. This might be a rea- 
sonable assumption for the potential of remote charged 
impurities located at dimp = 1 nm to 2 nm away, but not 
for vacancies or adsorbates that would be strongly local- 
ized on one of the layers. Without a microscopic theory 
that would give the matrix structure of the impurity po- 
tential in this chiral basis (and we note that such a theory 
would likely be non-universal depending strongly on the 
type of defect), it is more reasonable to set i^((/)) = 1 for 
short-range impurities. For generic short-ranged disor- 
der, this would correctly weight the relative importance 
of the interband and intraband contributions at the ex- 
pense of losing the chirality enhancement factor of order 
unity. Since we are interested in quantities such as the 
power-law trends for cr{n) (left panels of Fig.[2]and Fig.|3]) 
or the ratio between cr(r7,) at high and low carrier density 
(right panels of Fig. [2] and Fig. [S]), this approximation is 
well suited to the scope of this work. 



Appendix D: Bilayer graphene: Analytical results 

In this section we derive analytic results for the trans- 
port properties of bilayer graphene. The bilayer graphene 
Hamiltonian is 



H 



( vq-k* 

UOTT tj^ 

t± WoTT* 

V wqtt 



(Dl) 



where tt = ti(kx + ^^y) and energy eigenvalues and eigen- 
vectors (up to normalization) are 
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(D2) 



(D3) 



where 



+ v/(i±/2)2 + {hvoky 



(D4) 



From Eq. [5] taking into account only the low energy 
band with the energy e^, the Fermi velocity is given by 



matrix element into a plane-wave-like overlap Vq^^ and a 



_ 1 de 
^ h dk 



hvokp 

= Vn , 



(D5) 
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FIG. 5: (Color online) The chirality factor F{(f>) = j(s, k, <?!) = 
0|s,k',<^)P is determined by the wavefunction overlap be- 
tween initial and final states. Top panel shows ABA stacked 
trilayer graphene, and bottom panel shows ABCB stacked 
tetralayer graphene for Ef =0.1 eV. The bands are labeled 
according to the J-chiral decomposition shown in Table [H 
For intraband scattering, the chirality factor at small Fermi 
energy is given by Eq. IB3I However, the interband chirality 
factor is identically zero for scattering between the J — 1 and 
J — 2 bands of the ABA trilayer and strongly suppressed for 
scattering between the J — 1 and J = 3 bands of the ABCB 
tetralayer. As discussed in the text, we believe that this sup- 
pression is unphysical for generic short-range impurities. The 
insets in both panels show the tight-binding band structure 
determined numerically by solving Eq. [1] as discussed in the 
text and dashed lines indicate the Fermi energy. 



and the density of states per spin and valley at the Fermi 
energy is 



p(eF) 



2T:hvF 27r(?iuo)^ 



(D6) 



The chirality factor within the same band is given 

F(</>) = |(e+,0 = O|£+,0)|2 

- ^[l-V+il + v)cos4>]\ (D7) 



where 77 = l/^/T+n/no, n = gsg^k"^ / {'^tt) , m = 
9s9vko/{ATT) and hvoko = t±/2. 

For simplicity, let's consider the conductivity when the 
Fermi energy crosses only the low energy band e^. For 
short-range disorder, from Eq. [7] 



(D8) 



where f{ri) = |(5?7^ — 2ry + 1) with the chirality factor of 
Eq.[D3 while f{rj) = 1 for F{(t>) = 1. Then, we find 



h 



^imp 



/ hvo 



Note that if we include the inter valley scattering, the 
conductivity becomes smaller by a factor of 2. 
Similarly for the bare Coulomb disorder with dimp = 
we have 



a ■ 



1 



8rf 



3?7^ -277 + 3 



(DIO) 



while for screened Coulomb disorder with rfimp = we 
find (here a ^ 1) 

^2 



''imp 



5772 - 27/ + 1 



(Dll) 



These results are consistent with the numerical data 
shown in Fig. [5] and Fig. |3] as well as with earlier 
results 
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